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Abstract: It is anticipated that hard double parton scatterings will occur frequently in 
the collisions of the LHC, producing interesting signals and significant backgrounds to cer- 
tain single scattering processes. For double scattering processes in which the same hard 
scale t = ln(Q 2 ) is involved in both collisions, we require the double parton distributions 
(dPDFs) D 3 ^ 02 (xi, X2] t) in order to make theoretical predictions of their rates and prop- 
erties. We describe the development of a new set of leading order dPDFs that represents 
an improvement on approaches used previously. First, we derive momentum and number 
sum rules that the dPDFs must satisfy. The fact that these must be obeyed at any scale is 
used to construct improved dPDFs at the input scale Qq, for a particular choice of input 
scale (Ql = 1 GeV 2 ) and corresponding single PDFs (the MSTW2008LO set). We then 
describe a novel program which uses a direct x— space method to numerically integrate the 
LO DGLAP equation for the dPDFs, and which may be used to evolve the input dPDFs 
to any other scale. This program has been used along with the improved input dPDFs 
to produce a set of publicly available dPDF grids covering the ranges 
1(T 6 < x 2 < 1, and 1< Q 2 < 10 9 GeV 2 . 
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1. Introduction 

In the standard framework for calculating inclusive cross sections for hard scattering pro- 
cesses in hadron-hadron collisions, it is assumed that only one hard interaction occurs 
per collision (plus multiple soft interactions). This assumption is typically justified on 
the grounds that the probability of a hard parton-parton interaction in a collision is very 
small. Thus the probability of having two or more hard interactions in a collision is highly 
suppressed with respect to the single interaction probability. 

Hadron-hadron collisions in which two (or more) distinct pairs of partons hard scatter 
are nevertheless possible. Early theoretical studies of double scattering were carried out 
in the context of the parton model [1-3] , with subsequent extension to perturbative QCD 
[4-16] . Such processes have in fact been observed experimentally - both in y/s = 63 GeV 
pp collisions by the AFS collaboration at the CERN ISR [17] and more recently in y/s = 
1.8 TeV pp collisions by the CDF collaboration [18] and y/s = 1.96 TeV pp collisions by 
the DO collaboration [19] at the Fermilab Tevatron. 

The much greater energy and higher luminosity of the LHC implies that we will observe 
a greater rate of events containing multiple hard interactions in this experiment than in 
either of the two mentioned above. Moreover, a number of calculations [20-23] suggest that 
the products from multiple interactions will represent an important background to signals 
from the Higgs and other interesting processes. Further calculations [24-26] indicate that 
certain types of multiple interactions will have distinctive signatures at the LHC, facilitating 
a detailed study of this process by the experiment. 

The importance of multiple scattering signals and backgrounds at the LHC necessitates 
a good quantitative understanding of these processes. In particular, it is important to 
understand double scattering, which will be the dominant multiple scattering mode at the 
LHC. Assuming only the factorisation of the two hard subprocesses A and B, the cross 
section for this process in proton-proton scattering may be written as: 



The as are the parton-level subprocess cross sections. These are also encountered 
in single scattering, and are known for essentially all processes of phenomenological in- 
terest. The quantity m is a symmetry factor that equals 1 if A = B and 2 otherwise. 
The Tij(x\, X2,b;ti,t2) represent generalised double distributions. They may be loosely 
interpreted as the inclusive probability distributions to find a parton i with longitudinal 
momentum fraction x\ at scale t\ = ln(Qf) in the proton, in addition to a parton j with 
longitudinal momentum fraction X2 at scale £2 = h^Ql)) with the two partons separated 
by a transverse distance b. The scale t\ is given by the characteristic scale of subprocess 
A, whilst t2 is equal to the characteristic scale of subprocess B. Note that CDF and DO 
have measured double scattering via the 7 + 3jet final state, with A corresponding to 7+jet 
production and B to dijet production. 




i y j,k,l 



(1.1) 



x Tfcz (x\ , x' 2 , b; ti,t2)dx\dx2dx' l dx' 2 d 2 b 
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It is typically assumed that Fij(xi,x 2 ,b;t±,t 2 ) may be decomposed in terms of longi- 
tudinal and transverse components as follows: 

T ij (x 1 ,x 2 ,b;t 1 ,t 2 ) = Dii(x 1 ,x 2 ;t 1 ,t 2 )Fi(b) (1.2) 

The function D l £(x\, x 2 ; t±, t 2 ) has a rigorous interpretation in leading order (LO) 
perturbative QCD as the inclusive probability of finding a parton i with momentum fraction 
x\ at scale t\ and a parton j with momentum fraction x 2 at scale t 2 in the proton. Accurate 
prediction of double parton scattering cross sections and event signatures requires good 
modelling of D^(xi, x 2 ; t±, t 2 ) and Fj(b). In particular, one must correctly take account of 
the effects of correlations in both longitudinal momenta and transverse positions in these 
functions. 

Correlations between the partons in transverse space are highly significant - at the very 
least, they must tie the two partons together within the same hadron. As one might suspect, 
their precise calculation is not possible using perturbation theory. Existing models typically 
use Gaussian or exponential forms to describe the Fj(b), or a sum of Gaussian/exponential 
terms [27,28]. 

On the other hand, correlations in longitudinal momenta are typically ignored. The 
usual assumption (applied in the phenomenological calculations of [20-26]) is that at 
least for small Xi values the longitudinal momenta correlations are small, and therefore 
DV(xi,x 2 ;ti,t 2 ) may be taken to be equal to a product of the relevant single parton 
distribution functions (sPDFs) - i.e. DV(xi,x 2 ;ti,t 2 ) = D l h (xi;ti)D{(x 2 ;t 2 ). With this 
assumption, plus the assumption that F l -{b) is the same for all parton pairs ij involved 
in the double scattering process of interest, the cross section cr® AB ^ has the particularly 
simple form: 

^ = ~2^^~ (L3) 



0"cff 



d 2 b{F(b)f 



The quantity &( X ) ls ^he single scattering cross section for hard process X. The factor 
cr c ff in the denominator has the dimensions of a cross section. It can be understood as 
follows. Given that one hard scattering occurs, the probability of the other hard scattering 
is proportional to the flux of accompanying partons; these are confined to the colliding 
protons, and therefore their flux should be inversely proportional to the area (cross section) 
of a proton. Interestingly, the CDF and DO measurements give cr e g ~ 15 mb, which is 
roughly 20% of the total (elastic + inelastic) pp cross section at the Tevatron collider 
energy. 

It is argued that the approximation of the D l £(xi, x 2 ; t\, t 2 ) as a product of single 
PDFs should be particularly applicable in collider experiments where small x values are 
probed (i.e. large total system centre of mass energy with respect to subprocess energy) 
due to the large population of partons at these x values. The CDF experimental data also 
agree with this assumption, with no sign of x-dependence in their measured a e g over the 
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x ranges accessible to them (0.01 — 0.40 for their first subprocess, and 0.002 — 0.20 for the 
other) [29]. The DO data confirm this result, with no measured variation of a e s with the 
(second highest) jet transverse momentum. 

Even if the factorisation assumption holds at the CDF scale (i.e. Q 2 ~ 100—1000 GeV 2 ), 
it is unlikely that it will hold at higher scales such as will be encountered at the LHC. 
In [30,31], the behaviour of the distributions D l ?(x\,x 2 \t) = DV(xi,X2',t,t) with the two 
hard scales set equal (hereafter known as the dPDFs) were investigated. The authors de- 
rived an equation dictating the scaling violations (i.e. t dependence) of the dPDFs. This 
equation is an analogue of the DGLAP equation for sPDFs (sDGLAP equation) [32-35]. 
An important prediction of this equation is that, even if the dPDFs factorise at some 
scale to, then at any different scale factorisation will be violated [36]. In other words, 
the naive product D l h {x\; t)Di{x2\ t) where the D l h s satisfy sDGLAP is not a solution 
of the dDGLAP equations. Explicit numerical solutions of the LO 'double DGLAP' 
(dDGLAP) equation based on factorised inputs at Qq ~ 1 GeV 2 suggest that the vi- 
olations are significant even for small x, with deviations on the order of 10 — 30% at 
Xl = x 2 ~ 0.1, Q 2 ~ 10 4 GeV 2 [37,38]. 

On a more fundamental level, the factorised approach is inadequate in that it fails 
to take account of even very basic correlations associated with the fact that finding a 
quark of given flavour reduces the chances of finding another with the same flavour. It also 
neglects correlations associated with the fact that finding a parton with momentum fraction 
x reduces the probability of finding another parton with momentum fraction near 1 — x 
(apart from a crude cut-off in the probability when the sum of the two parton momenta 
exceeds 1). In other words, factorised forms do not obey the relevant momentum and 
number sum rules (see section 3). In [27], a method to construct double (and multiple) 
parton distributions is suggested which is intended to capture some of the main features of 
these flavour and momentum correlations. However, this method does not introduce the 
correlations fully rigorously using the sum rules. Furthermore, the dPDFs are constructed 
entirely out of sPDFs, and no use is made of the dDGLAP equation in this paper. 

For the purposes of the accurate prediction of double parton scattering cross sections 
and signals at the LHC, there is a need for a more theoretically sound set of double 
distributions than either the naive factorised forms traditionally used, or the improved 
forms suggested in [27]. Here, we have attempted to address this issue for the specific 
case of the dPDFs (with t\ = t 2 ). First, we derive sum rules corresponding to momentum 
and valence quark number conservation that the dPDFs must satisfy. These are used as 
an aid to construct 'improved' dPDFs at the scale Qq = 1 GeV that correspond to the 
MSTW2008LO sPDF inputs [39]. The low scale dPDFs are then used as an input in a 
program we have written which numerically integrates the LO dDGLAP equation to higher 
scales. The end products of this paper are a set of LO dPDF grids covering the ranges 
10~ 6 < xi < 1, 10~ 6 < x 2 < 1, 1 < Q 2 < 10 9 GeV 2 , and all possibilities for the parton 
indices i and j. These grids, in addition to a simple interpolation subroutine designed to 
extract from the grids a dPDF value at a given x\,x 2 and Q, can be found at Ref. [40]. 

This paper is organised as follows. We begin with a brief review of the dDGLAP 
equation in Section 2. The dPDF sum rules are introduced and discussed in Section 3, where 
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we also explain how we have used these rules to construct input dPDFs at Qq = 1 GeV 
corresponding to the MSTW2008LO sPDF inputs. In Section 4, the details of the numerical 
procedure designed to evolve the input distributions to higher scales using the LO dDGLAP 
equation are given. Section 5 examines the ways in which our dPDFs differ from those 
obtained using previous approaches. Finally, we conclude in Section 6 with a summary 
and discussion of potential future directions for the work. 

2. The Double DGLAP Equation 

It is well established that in QCD, the parton content of the proton that is observed by 
a hard probe with virtuality Q 2 (or more than one hard probe with this virtuality) is 
dependent on the size of the virtuality. This dependence is explained by the fact that a 
harder probe (with a shorter associated wavelength) is able to 'see' finer scale structure 
in the proton, and in particular is able to resolve parton splittings that were unresolvable 
using a lower t = ln(Q 2 ) probe. This implies that parton distributions must be dependent 
on the scale t at which the proton is probed. There is a shift of these distributions towards 
smaller x values as t increases, as a consequence of a greater number of splittings being 
resolved. 

One can visualise the change in parton distributions as one probes at steadily higher 
scales t than the low scale to as a spacelike branching process originating from the initial 
distributions at scale to- As one probes to higher scales, one effectively progresses further 
in the branching process. 

The dDGLAP equation is a renormalisation group equation describing the change of 
the dPDFs with the hard scale t. It is based on the leading logarithm approximation (LLA) 
of perturbative QCD (the same is true for the sDGLAP equation). This approximation 
corresponds to a picture of the spacelike parton branching process from the low scale to to 
the probe scale t in which gluon emissions along the parton branches are strongly ordered 
in transverse momentum. Gluons emitted 'earlier' in the branching process are restricted 
to have smaller transverse momenta than those emitted closer to the probe scale. The 
dDGLAP equation effectively resums leading powers of [a s t] n generated by these gluon 
emissions to give the dPDFs at scale t. 

In [30,31], the following form for the dDGLAP equation is derived: 



Technically, the argument t in the above is the factorisation scale (which in practical 
calculations is typically set equal to the characteristic hard scale of the subprocesses) . The 



dDl 1 ^(x 1 ,x 2 ;t) 
dt 



a s {t) 
2vr 





(2.1) 
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renormalisation scale has been set equal to the factorisation scale to obtain this equation 
(as is conventional in a leading order analysis). 

In addition to the dPDFs and sPDFs D 3 h (x; t), the equation (2.1) contains two different 
types of splitting functions. The first are the well-known splitting functions P{^j(x) previ- 
ously encountered in the context of the sDGLAP equation. They are given to both LO and 
NLO in [41]. At leading order, the function Pi^j{x) may be interpreted as the probability 
of a parton i splitting to give a parton j with a fraction x of the longitudinal momentum 
of the parent parton and a transverse momentum squared much smaller than Q 2 (where 
t = ln(Q 2 )) [35]. The second, the Pi^jk(x), are new. They may be interpreted at LO as 
the probability of a parton i splitting to give the two partons j and k, the first of which 
has a fraction x of the linear momentum of the parent parton, the second of which has the 
remainder of the linear momentum 1 — x, and both of which have transverse momentum 
squared much less than Q 2 . 

The splitting functions Pj_>.j(x) each possess a large negative contribution at x = 1 
(these are contained within the 'plus prescription' functions together with explicit delta 
functions in the definitions). This contribution is included to take account of the fact that 
splittings of the parton i into other partons with lower momentum act to reduce the pop- 
ulation of partons with the original momentum. At a fundamental level, the contributions 
at x = 1 result from virtual gluon radiation diagrams. 

On the other hand, the functions Pj_^-fc(x) do not contain such contributions. This is 
to be expected as a virtual process is clearly not able to achieve the 1—^2 splitting i — > jk. 
At LO, the function Pj_^(x) is related to the 'real splitting' part 1 of the normal splitting 
functions P^-(x) according to: 

P^(x) = Y. P ^ x ) ( 2 - 2 ) 

k 

Equation (2.2) is the simple statement that the probability of splitting i — >■ j + anything 
is equal to the sum of probabilities of splitting i — > j + k, summed over all possibilities for 
k. 

A further simplification to (2.2) is possible at LO. Due to the fact that QCD only 
allows certain types of three particle vertices (i.e. triple gluon vertices and 'gluon emission 
from a quark' type vertices), the LO Pj_>jfc(x) is only nonzero for a small number of {i, j, k} 
combinations. In fact, given i and j, there exists at most one choice for k which makes 
Pi^jk{x) nonzero. We shall denote this special value of k by For example, 

is g when i = q, h j = and % when i = g,j = 

Given this fact, we note that (2.2) must contain at most only one term on the right 
hand side, and we may write: 

P&j(x) = I\^ Mi>j) (x) (2.3) 

In (2.3), we have extended the definition of to cases where there exists no choice 

for k to make P^j^x) nonzero. In these cases, can be chosen to be any parton, as 

both the right and left hand sides are zero for any choice. 

lr rhe functions P[^j(x) are obtained from the functions Pi^j(x) by dropping the terms proportional to 
6(1 — x). This includes removing plus prescription + signs where they appear. 
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Equation (2.3) effectively defines P^jk for all cases in which it is nonzero. At LO 
then, we may construct the following definition for P^jk- 

Pi , l (.)4^ ( " ""^ (2.4) 
[ otherwise 

It is interesting to consider the generalisation of the P^jk{x) functions to NLO (and 
indeed higher orders). Here one encounters a problem, in that the function P^jk(x) only 
has one longitudinal momentum argument because it is assumed that the parton k must 
possess the remaining longitudinal momentum originally carried by i that was not given 
to j, i.e. 1 — x. This is certainly true at leading order, where only two partons can be 
produced in a single splitting, by conservation of momentum. However, it is not true in 
general at NLO, where a single splitting can contain two QCD vertices, and produce three 
partons. At NLO and above, the splitting function P^jk should have two arguments, x\ 
and x 2 . 

The expansion of the more general function P^jk(x\, x 2 ) in terms of powers of a s 
would read as follows: 

P i ^ jk (x 1 ,x 2 ) = 5(1 - X! - x 2 )P ( >%{x 1 ) + ^pW jk ( Xl ,x 2 ) + . . . (2.5) 

The higher-order coefficients in this expansion cannot be obtained trivially from the higher- 
order coefficients of the splitting function Pi^j(x) as in the LO case. The general relation 
between the two for x < 1 is: 

= E f X1 dX2 p ^d x ^) ( 2 - 6 ) 
k J ° 

The unintegrated P^j k (x\,x 2 ) must contain more information than the P^j(x) for n > 0, 
and hence cannot be obtained from them. 

A consequence of the fact that Pi^,j k {x) with a single longitudinal momentum argu- 
ment x is not the right function to use at NLO and above is that the precise structure of 
the dDGLAP equation given in (2.1) can only be applicable at LO. 2 In what follows we 
restrict our analysis to the LO case, but will return to the generalisation of (2.1) to NLO 
in a future study. 

One can interpret the terms on the right-hand side of (2.1) using the parton branching 
picture. 3 Consider the inclusive probability of finding a pair of partons in the proton with 
flavours j\ and j 2 and longitudinal momentum fractions between x\ and x\+5xi and x 2 and 
x 2 + 5x 2 respectively at scale t, Dj^- 12 (xi,x 2 ;t)5xi5x 2 . It is obvious that when t is increased 
to t + At, two types of process may contribute to the change in this quantity. Splittings 
from higher-momentum partons giving rise to jij 2 pairs with the correct momentum act 
to increase this quantity, whilst splittings within the jij 2 pairs to give partons of lower 
momentum act to decrease the quantity. 

2 This is not explicitly stated in [30,31]. 

3 We use similar arguments as are used in Section 5.2 of [41] to explain the terms on the right hand side 
of the sDGLAP equation. 
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x 2 ;i)<5x 2 




Figure 1: Splitting processes that increase the population of ju'2 parton pairs with momenta 
in the ranges xi — ► xi + &ci, X2 — > x 2 + 5x2- P^j{x) is the 'real splitting' part of the 
splitting function P^Ax) - i.e. the splitting function minus the terms proportional to 
5(1 -x). 



At leading order in a s (which, as we have established, is the order under which equation 
(2.1) was derived), there are three types of splitting process that give rise to a pair of partons 
J1J2 with momenta in the ranges x\ — > x\ + 5x±, x 2 — > X2 + <5x 2 . The three are drawn 
schematically in Fig. 1. 

In the first, we start with a pair of partons j[j2 with momenta in the ranges x' x — > 
x[ + Sx'i, X2 — > X2 + 5x2- The quantity x[ must satisfy xi < x[ < 1 — x 2 - i.e. be large 
enough that j[ can split to give ji, and be small enough that the initial pair of partons is 
not carrying more momentum than the proton they are in. The parton j[ then splits to 
give as one of the products a ji with momentum in the range xi — > x\ + bx\. The second 
process is very similar but involves a splitting in the second parton. The third involves a 
single parton j' with just the right momentum x\ + x 2 — > x\ + x 2 + 6x2 splitting to give as 
its two daughters the pair jij 2 with momenta in the appropriate ranges. 

The leading order splitting processes reducing the population of partons with 
momenta in the given ranges are given in Fig. 2. There are two processes - in the first, the 
ji parton splits to give lower momentum partons, whilst in the second the j 2 splits. 

The correspondence between the diagrams of Figs. 1 and 2 representing the splitting 
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D 3 h 132 (xi,x 2 ;t) 5x x 5x 2 




D h (xi,x 2 ;t) 5xiSx 2 



+ 



r 




D 3 ^ 32 (xx,x 2 ;t) 5x^X2 



Figure 2: Splitting processes that decrease the population of j±j 2 parton pairs with mo- 
menta in the ranges x\ — > x\ + 5x\, x 2 — > x 2 + 5x 2 . Pj-^j is equal to the sum of the 
coefficients of the 5(1 — x) terms in the splitting function Pa^Ax) (including 5(1 — x) terms 
contained within plus prescription functions). 



processes and the terms on the RHS of (2.1) is fairly clear. Suitable labels have been 
added to the figures to bring out this correspondence. It is important to note that the first 
two sets of terms on the RHS of (2.1) cover four diagrams - the first two of both Fig. 1 
and Fig. 2. The 'real splitting' parts of the terms correspond to the diagrams in Fig. 1, 
whilst the 'virtual correction' parts correspond to the diagrams in Fig. 2. We also note 
that the last set of terms on the RHS of (2.1) contains sPDFs because it corresponds to 
the diagram in which a single parton splits to give the pair j\j 2 . There is no integral in 
these terms because of the property of LO QCD that a single splitting can only give rise to 
two partons. Thus the single parton that splits is essentially restricted to have momentum 
exactly equal to x\ + x 2 . We shall hereafter refer to the last set of terms on the RHS of 
(2.1) as the 'sPDF feed' terms, for obvious reasons. 

A solution to (2.1) in terms of sPDFs is obtained in [30,31], and presented for the first 
time in x-space in [36]. Let us introduce the 'natural' evolution variable r defined in terms 
of t according to: 



to 
1 

2^b 



In 



2tt 
" t - 



(2.7) 



ln(A 2 
to ~ ln(A2 



QCD) 



at LO 



QCD J 
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In terms of the variable r, the solution to the dDGLAP equation reads: 

Di^(x 1 ,x 2 ;r) = Z^ c 2 orr) (*i, x 2 ; r) (2.8) 

f 1 -* 2 dzi ( x ~ zx dz 2 n j'j> 

where: 

. . FT I-1—X2 A-,^ fI—zi J 

K&rr^^r) = Y, / W |i / ( Zl + Z2 ; T>) (2.9) 

Jo J xx z l Jx 2 z 2 

J JXJ2 

Z!+Z 2 \Z 1 +Z 2 J * / / 

The Green's functions (x; r) are defined such that they satisfy the initial conditions 
Dj (x; r = 0) = <5jj5(l — x) and change with r according to the sDGLAP equation: 

In effect, the function Z}? (x; r) gives the inclusive probability that one finds a parton j 
with longitudinal momentum fraction x at scale r inside a dressed object that looks like a 
pure i parton at the scale r = 0. 

A pictorial representation of the solution (2.8) in terms of parton branching is given 
in Fig. 3. One observes the need to specify some initial conditions D 3 ^ J2 (xi,x 2 ;t = 0) to 
obtain the distributions at higher scale, which is a direct reflection of the fact that the 
dDGLAP equation can only predict changes in the distributions with r. 

The depiction of dDGLAP evolution as in Fig. 3 leads us to make a suggestion as to 
how one might calculate the double distributions for which the two scales are not equal, 
D % £ (xi, x 2 ; ti,t 2 ). The arguments T\ and t 2 in this distribution correspond to the factorisa- 
tion scales for parton i and j respectively. Consider the analogous figure to Fig. 3 for these 
distributions. It seems likely that this figure would be the same, except with t\ replacing r 
on the 'upper legs' of the diagrams, t 2 replacing r on the 'lower legs' of the diagrams, and 
the upper limit of the r' integration replaced by min(ri,T2). If this ansatz is correct, the 
double distributions D^(xi, x 2 ; t\, t 2 ) with (say) n < t 2 should be calculated by taking 
the dPDFs with r = t±, and then performing sDGLAP evolution at each x\ from t\ to t 2 
in the x 2 variable. The upper limit in the sDGLAP evolution at given x\ should be 1 — x\. 
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Figure 3: A schematic representation of the solution of the dDGLAP equation (2.8) 
terms of the parton branching picture. 
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3. The Double Parton Sum Rules and the Initial Distributions 
3.1 The Double Parton Sum Rules 

It is well known that the sPDFs satisfy two types of sum rules which represent the fact 
that both momentum and valence quark number should be conserved under evolution. One 
might wonder whether corresponding rules exist for the dPDFs. We have managed to prove 
that the following equalities are preserved by dDGLAP evolution, provided they hold at 
the starting scale io : 



Momentum Sum Rule: 

Let M be the momentum fraction carried by the proton (= 1). Then: 

V f 1 ^ dx l x 1 D{ lj2 (x l ,x 2 ;t) = (M - x 2 )D{ 2 (x 2 ;t) (3.1) 

Number Sum Rule: 

Let ji v = j\ — ji (ji 7^ g), and Nj lv be the number of 'valence' j\ quarks in the proton. 
Then: 



rl—X2 

/ dx 1 D 3 ™ 2 (x 1 ,x 2 ;t) = < 
Jo 



N hv D h ( x 2 ; t ) when j 2 / j x or j 1 

(N jlv -l)Di*(x 2 ;t) whenj 2 =n (3.2) 
[(N jlv + l)Df(x 2 -t) wheni 2 = J 1 



The only nontrivial inputs to this proof are the following relations, which must be 
obeyed by the splitting functions in order that the number and momentum integrals are 
conserved for the sPDFs: 

V" / dxixiPj'^j (xi) = 0; / dxiPj'^ jv (xi) = (3.3) 
j Jo Jo 

We do not present the proof here, since it is straightforward and rather lengthy. A crucial 
point we must emphasise however is that the precise structure of the sum rules given above 
is required in order that they should be preserved by dDGLAP evolution. In particular, 
the prefactors in front of the sPDF quantities on the right hand sides of (3.1) and (3.2) 
must have the values given. 

By appropriately combining equations (3.1) and (3.2) with the sPDF momentum and 
number sum rules, one can construct integrals over both arguments of the dPDFs which 
give conserved quantities such as M or Nj v (or products of these quantities). Examples of 
such integrals are given below: 



V/ dx 2 dx 11 ^-D J h lja (x 1 ,x 2 ;t) = M = l (3.4) 

^ Jo Jo M — x 2 



L 



J1J2 
I fl-X-2 



Dr(xi,x 2 ;t) D{^(x 1 ,x 2 ;t) 
'o^Jo ^^J—l N-T^ = N ^ (3 ' 5) 
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These relations are preserved under dDGLAP evolution. By contrast, integrals such 



as X^ij 2 lo dx i Jo X2 dxiXix 2 D{ lj2 (xi,x 2 ;t) and dx 2 Jq X2 dxiD 3 h lvjlv (xi,x 2 ;t), which 



one might naively think should give conserved momenta or valence quark numbers, are not 
conserved by dDGLAP evolution and so do not correspond to such physical quantities. 

An appealing interpretation of (3.1) and (3.2) exists in terms of probability theory 
(although such a picture has in no way been used to obtain these relations). The dPDF 
sum rules £1X6 clllclls gous to the result in probability theory that for two continuous random 
variables X and Y, the probability density functions relating to X and Y must satisfy 4 : 



The integral is performed over all values that X can take given that Y = y, and 
E(X a | Y = y) is the expectation value of X a given that Y has value y. All of the prefactors 
on the right hand sides of Eqns. (3.1) and (3.2) are essentially conditional expectations as 
in (3.6). The (1 — x 2 ) factor on the right hand side of (3.1) is the conditional expectation 
value for the momentum of all of the other partons in the proton given that one has found 
a parton of longitudinal momentum fraction x 2 . The (Nj lv — 1) factor for the j 2 = j\ 
case of (3.2) is the conditional expectation for the number of j\ partons minus the number 
of ji partons elsewhere in the proton, given that one has found a parton of flavour j±. 
The prefactors for the other number sum rule cases may be interpreted as conditional 
expectation values using similar logic. 

The fact that the forms of the number and momentum sum rules can be justified 
using general arguments strongly suggests that these rules should hold to all orders in 
perturbation theory, just as the sPDF number and momentum sum rules do. We have 
been restricted to an LO proof that they hold at all scales if they hold at the starting scale 
by the fact that we only have the LO dDGLAP equation. 

Although we have not derived the momentum and number sum rules from first princi- 
ples, they appear to satisfy a number of non-trivial consistency checks. For example, one 
might worry that there might not be a set of dPDFs that satisfy the full set of rules. A 
potential source of tension between the different rules is the integral: 



One can evaluate this integral using (3.1) or (3.2), in combination with appropriate 
sPDF sum rules. If different results were produced depending on which dPDF sum rule 
was used, this would indicate an inconsistency. However, one obtains the same result, 
Nj lv — fj lv (where fj lv is the momentum fraction carried by valence j\ partons), with 
either approach. This lends further credibility to the sum rules (3.1) and (3.2). 

4 One should bear in mind that the correspondence between (3.6) and (3.1)/(3.2) is not completely 
straightforward, as the parton density functions are not really simple probabilities. Rather, they may be 
better interpreted as number distributions. This results in, for example, the sPDFs being normalised to 
the number of partons of the given type in the proton rather than 1. Of course an analogous relation to 
(3.6) exists for such distributions. 




(3.6) 




(3.7) 
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It is notable that the complete set of dPDF sum rules, (3.1) and (3.2), do not appear 
anywhere in the extant literature, although similar sum rules have been derived for the two- 
particle fragmentation functions in [42]. An early paper on the subject, [3] (see also [10]), 
introduces some 'constraints' resembling the number sum rules, which are used as an aid 
in constructing some simple model dPDFs. However, the constraints are only imposed for 
two specific dPDF cases, and the paper does not make any explicit statement about the 
general form of the number sum rule. In particular, they do not describe the subtleties of 
the number sum rule with regard to the different possible proportionality constants on the 
right hand side of (3.2). 

In some sense, the dPDF sum rules are more restrictive than their sPDF counter- 
parts. The sPDF sum rules state that the quantities M = ^ L 1 dxxD l h (x; t) and Ni v = 
Jq 1 dxD^(x; t) are conserved under evolution whatever their initial values, and we make the 
physical choices M = l,N Uv = 2,Nd v = 1 for the proton. On the other hand, Eqns. (3.1) 
and (3.2) are only preserved under evolution if they hold at the starting scale. This is 
linked to the fact that one initially has the freedom in the sum rules to specify the momen- 
tum/parton composition of the hadron M and N iv (although M ^ 1 is not very physical). 
However, once these have been specified in the sPDF sector, the structure of the multipar- 
ton sum rules is effectively fixed. 

The restrictive nature of the dPDF sum rules can be used to place nontrivial constraints 
on the input distributions that are physically allowable in the dDGLAP equation. If we 
believe that the dPDF sum rules should hold at the starting scale, then we can use the 
constraints provided by the rules to improve on the factorised inputs previously used at 
the starting scale Qq ~ 1 GeV 2 . This is discussed in the next section. 

3.2 Use of the Double Parton Sum Rules to improve the Input Distributions 

As was mentioned in Section 1, it is a common assumption that the input double distri- 
butions should be equal to the product of the relevant sPDFs at low X\ and x^- The logic 
behind this is that there exist large populations of partons of all active flavour types and x 
values at low x. Given these large populations, we would expect the extraction of a parton 
with a given flavour type ji and small longitudinal momentum x\ not to have a strong 
effect on the probability of finding another parton of flavour ji (where j2 can be equal to 
ji) and small longitudinal momentum X2- This leads to a joint probability distribution 
which can be expressed as a product of single distributions at low x\, X2- 

This factorisation assumption appears to be backed up by the available CDF and DO 
data. Consequently, we would like our improved input dPDFs to maintain a factorised 
form for low X\,X2, whilst now obeying the sum rules (3.1) and (3.2). The first question 
to be addressed in this section is whether this is in fact possible for all the dPDFs, i.e. 
whether the sum rules are compatible with factorisation at low x\,X2 in all cases. 

To help answer this question, we introduce the 'double evolution' representation for the 
dPDFs. In this representation, the well-known {singlet, gluon,valence,tensor} / {£, g, Vi,Ti} 
combinations (defined in, for example, Chapter 4 of [41]) are used as the flavour basis for 
both parton indices in the dPDF. The relationship between this basis and the 'double 
human' basis in which both parton indices i,j are one of g,u,u etc. can be clarified using 
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dPDF Type 


Relevant Sum Rules 


Valence- Valence 
Valence- Sum 
Valence- Tensor 
Tensor- Tensor 
Tensor-Sum 
Sum-Sum 


Number (involved in two rules) 
Number + Momentum 
Number 
None 

Momentum 

Momentum (involved in two rules) 



Table 1: The different dPDF classes under the 'double evolution' representation of the 
dPDFs, and the types of sum rules each is engaged in. 



an example: 



jyT 3 u v _ jj{u+u-d-d){u-u) 
h h 



(3.8) 



ryuu 



r\du 
U h 



7-1 du 
u h 



7-1 UU 
U h 



du 



The longitudinal momentum arguments of each term in this equation are the same. 
The use of the 'double evolution' representation has the advantage that it splits the dPDFs 
into six sets, each of which must satisfy different combinations of the sum rules. We refer 
to the singlet and gluon combinations as the 'sum' combinations (as they describe the sum 
of quark and gluon contributions respectively). Since ^ j = £ + g, any dPDF with a 
'sum' flavour index will be involved in a momentum sum rule, whilst any dPDF with a 
'valence' flavour index will be involved in a number sum rule. Those dPDFs where each 
of the indices are one out of the 'sum' and 'valence' combinations will be involved in two 
sum rules. 

The six sets of dPDFs along with the combinations of sum rules each is involved in 
are given in Table 1. We do not write out the explicit forms of the sum rules under the 
double evolution basis in this table. To obtain each rule, one must first construct the 
appropriate integral (i.e. J dxiXi[D^ k (xi,x 2 ) + Df i k (xi,x 2 )] for a momentum sum rule or 
J dx\X\D % £ k (x\, x 2 ) for a number sum rule, where k can be any double evolution basis 
index). The sum rule is then obtained by expanding each dPDF in the integral in terms 
of human basis dPDFs (as in (3.8)), followed by the use of equations (3.1) and (3.2). We 
illustrate this procedure for the case of the u v T% number sum rule: 



rl-xi 

Jo 



dx 2 D^( Xl ,x 2 ) 



1 — x± 



dxo 



D^( Xl ,x 2 ) + Dl^( Xl ,x 2 ) 



(3.9) 



D^( Xl ,x 2 ) 



D d ^{x u x 2 ) 



={N Uv - l)D u h ( Xl ) + {N Uv + l)DUxi) 

- N Uv D d h ( Xl ) - N Uv D d h ( Xl ) 
=N Uv D^( Xl )-D^( Xl ) 

If one investigates each class of dPDF and their respective sum rules, then one finds 
that in most cases one is allowed dPDFs which satisfy the sum rules and are approximately 
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equal to the product of single distributions at low x\ and x 2 . There is however a type of 
dPDF for which these two requirements cannot be simultaneously satisfied - the dPDF 
with two of the same valence combinations as its flavour indices (e.g. D% vUv ). 
The number sum rule that this type of dPDF must satisfy reads: 

pl—X2 

/ dx 1 D{^{x 1 ,x 2 ;t ) = N jv D%{x 2 ;t ) - D{ +J {x 2 ;t ) (3.10) 
Jo 

Consider this equation for small x 2 . Assuming no pathological behaviour of the function 
Dj^ Jv (xi, x 2 ; to) near the kinematical bound X\ + x 2 = 1, the integral on the left hand side 
of (3.10) is dominated by contributions from the small x\ region where D^ 3v (xi, x 2 ; to) 
is largest. A factorised form for D : t Jv (xi, x 2 ; to) at small x\,x 2 would then result in the 
left hand side behaving like x^ (where x~ av is the small x behaviour of a typical valence 
sPDF). 

On the other hand, the right hand side of (3.10) is dominated by the — Dj^~ J (x 2 ; to) 
term. This is due to the fact that this term receives contributions from the sea, and sea 
sPDFs diverge faster than valence sPDFs at low x. We expect —Dj^~ J (x 2 ) to behave like 
— x^" (where a typical sea sPDF behaves like x a " at low x). The right hand side then 
behaves very differently 5 from the left hand side, and it is impossible to satisfy the sum 
rule (3.10) using a dPDF that factorises at low xi,x 2 . 

We conclude that we must abandon the possibility of factorisation into a product of 
sPDFs at low x\,x 2 for the D^ 3v {x\, x 2 ; to). The fundamental origin of the second term 
on the right hand side of (3.10) which precludes the possibility of a factorised form for 
Dj^ 3 " (x\, x 2 ;to) is of course in number effects. By 'number effects' we mean the fact that 
finding a parton of a given type alters the probability of finding a further parton of the 
same type, due to the fact that the number of that parton has decreased. 

The CDF and DO results are not in contradiction with the above conclusion, since in 
these experiments the vast majority of double parton scatterings observed would have been 
initiated by gluons and sea quarks. The dPDFs relevant to these partons are able to have 
factorised forms at low x\,x 2 . 

At first glance, it might appear that the statement of the inadequacy of factorised forms 
as applied to the valence- valence distributions has already been made, in [10]. However, 
our statement and the one in [10] are really very different things. In [10], the authors 
argue that one should not use a factorised form for the valence-valence dPDFs at large 
xi,x 2 . The reasoning behind this is that the inaccuracies of the factorised ansatz at large 
x±,x 2 due to the fact that it neglects momentum conservation effects are most strongly 
noticed in the valence- valence dPDFs, which are dominant at large x\, x 2 . Whilst we agree 
with their conclusions, we further propose that the factorised forms should not be used to 
describe equal flavour valence- valence dPDFs at small x\,x 2 , a point that is missed in [10] 
and elsewhere. 

Bearing in mind the points made above, we proceed to discuss how some input distribu- 
tions approximately obeying the sum rules might be obtained. One might initially wonder 

5 Regge theory arguments, for example, would suggest a v ~ \ and a s ~ 1, and 'modern' global fit sPDFs 
show a similar trend. 
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whether it is possible to develop a framework for constructing dPDFs out of combinations 
of sPDFs that does not make reference to any specific choices for the input sPDFs (e.g. 
MSTW, CTEQ). Instead, it would make intelligent use of the sum rules the sPDFs have 
to satisfy to ensure the dPDF sum rules were satisfied. However, investigation into this 
route has revealed that a framework of this kind does not seem to exist, even to construct 
dPDFs that only satisfy one of the two types of sum rules. 

Our discussion must therefore be based around some specific set of input sPDFs. For 
the purposes of producing the most accurate set of dPDFs we can, it would seem sensible 
to use the inputs from the most recent LO fit by one of the PDF fitting collaborations. 
We have chosen to use a set which almost exactly corresponds to the MSTW2008 LO 
inputs (Equations 6-12 and the first column of Table 4 in [39], with Qq = 1 GeV and 
a s(Qo) = 0.68183). The only differences between our inputs and those of [39] are that we 
have set the initial s v distribution to zero, and have added the following terms to the d 
distribution: 

-148.103388x 3 (l - x) 10 ' 8801 + 500x 4 (l - x) ia8801 (3.11) 

These modifications have been made in order to fix the problem that the MSTW2008 
LO s and d input distributions go slightly negative in some region of x. Even though strictly 
speaking these LO sPDFs should never go negative, the deviations below zero observed in 
the MSTW2008 LO s and d inputs are perhaps tolerable in single scattering calculations 
due to their small size (s, d > —0.0005). However, we must insist on using sPDFs which 
are strictly non-negative when expressed in the 'human' flavour basis 6 to build our input 
dPDFs. We can explain why this has to be the case by considering the dPDFs in the 
'double human' basis in which at least one flavour index corresponds to an sPDF which 
goes negative. Like all LO dPDFs in the 'double human' basis, they cannot go negative 
(due to their interpretation as a probability) . If we use a pseudo-factorised prescription to 
construct the dPDFs, then these dPDFs will go very seriously negative where the sPDF 
in one direction takes small negative values, and the sPDF in the other becomes large and 
positive. We therefore require strictly non-negative input sPDFs. 

We can identify two key features that we would like to build in to our set of input 
dPDFs. These are the following: 

1. The dPDFs should be suppressed below factorised values near the kinematical bound 
(i.e. the line X\ + %2 = 1) due to phase space considerations. 

2. Terms should be added/subtracted from certain dPDFs to take account of number 
effects. 

Let us begin by discussing how the first requirement might be incorporated. In the 
early papers [3, 5, 7, 10], a common (1 — x\ — £2) suppression factor multiplying all of 
the dPDFs was advocated. This was motivated by arguments based on the recombination 
model of [43], or the Kuti-Weisskopf model of [44]. More recently [38], it has been suggested 
that a higher power of (1 — x\ — X2), such as (1 — x\ — X2) 2 , might be appropriate. With 

6 The 'human' flavour basis is the one in which the parton index i = g,u, u, etc. 
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the benefit of knowledge of the sum rules, we can see that neither of these alternatives is 
entirely satisfactory. To illustrate this, let us just consider the momentum sum rule for 
the moment (which is the relevant rule with regards to phase space considerations), and 
let us consider the lines x\ = and X2 = 0. Along these lines, all momentum sum rules 
are perfectly satisfied using factorised dPDFs, whilst dPDFs including a (1 — x\ — x-i) or 
(1 — x\ — X2) 2 factor violate the sum rules badly. 

Thus a (1— x\ —X2) n factor alone multiplying all of the dPDFs suppresses the functions 
rather too severely near the lines x\ = and X2 = 0, and it would seem that a phase space 
factor which approached 1 near these lines would be more desirable. We can actually make 
sense of this from an intuitive point of view. The phase space suppression factor is inserted 
to take account of the fact that finding a parton with x = x\ reduces the probability of 
finding another parton with x = X2 if x\ + X2 is close to 1. One would expect a much 
smaller reduction if x\ were small and X2 were large than if both x\ and X2 were large, 
even if the sum of x\ and X2 was the same in both cases. Indeed, one would anticipate that 
the reduction should tend to zero as X\ (or X2) tended to zero - that is, the phase space 
factor should approach 1 as one approaches the lines x\ = and X2 = 0. 

Here, we continue to follow the tradition set by previous papers in that we have at- 
tempted to apply a universal phase space factor to all of the dPDFs. Use of a (positive) 
universal phase space factor has the advantage that it is guaranteed to produce positive 
double human basis dPDFs. However, instead of using (1 — x\ — X2) n alone, we tried the 
following as a 'first guess' for the phase space factor p, motivated by the above discussion: 



Following the more recent work by Korotkikh and Snigirev [38], we choose n to be 2. 
This choice of phase space factor gave dPDFs which satisfied the momentum sum rules 
reasonably well. In the left panel of Fig. 4, we plot the 'sum rule ratio' with this phase 
space factor for the particular example of the (X + g)g momentum sum rule - the sum 
rule ratios for the other momentum sum rules exhibit very similar behaviour. The sum 
rule ratio for a particular sum rule and set of dPDFs is defined as the sum rule integral 
calculated using the dPDFs divided by the sPDF quantity it should be equal to. It is a 
function of an x variable, and measures how well the dPDFs satisfy the given sum rule - 
the closer the ratio is to 1 over the full x range, the better the dPDFs satisfy the sum rule. 

On the other hand, the dPDF number sum rules are not particularly well satisfied by 
this prescription (this is illustrated in the right panel of Fig. 4). This is true even for those 
dPDFs which are involved in a number sum rule but which are not affected by number 
effects - e.g. u v d v . For these dPDFs, the phase space factor alone should be sufficient 
to cause the dPDFs to satisfy the relevant number sum rules - thus our first guess is not 
fully satisfactory. We have discovered that a slight adjustment to the form (3.12) resolves 
this problem. Let us allow the phase space factor to depend on the parton indices i,j on 
the dPDF such that (prior to adjustments relating to point 2 above) the input dPDFs are 
constructed according to: 



P(X 1 ,X 2 ) = (1 - X! - X 2 ) n {l ~ Xl) "(1 - X 2 ) 



—n 



(3.12) 



D%{x 1 ,x 2 ;t ) = D l h (x 1 ;t )Dl(x2;t )p l3 (xi,x 2 ) 



(3.13) 
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aT 0.6 



Figure 4: Sum rule ratios for the (S + g)g momentum and u v d v (integrating over u v ) 
number sum rules, when the phase space factor is as given in (3.12) with n = 2. 



We now define p l i{x\,Xi) as follows: 

p ij ( Xl ,X 2 ) = (1 - X! - X 2 ) 2 (l - Xl)- 2 - Q ^(l - X 2 y 2 - a{i) 

where: 



a(i) 



if i is a sea parton 
0.5 if i is a valence parton 



(3.14) 



(3.15) 



If either i and/or j contain both valence and sea contributions, then one should con- 
struct the dPDF by taking the factorised product, splitting it into sets of terms corre- 
sponding to valence- valence, valence-sea, sea-sea, etc., and then applying the appropriate 
phase space factor to each set of terms. Note that the phase space factor is no longer 
universal, but is nearly so - it turns out that this prescription is guaranteed to produce 
positive human basis dPDFs provided all the valence sPDFs are positive, which is the case 
for the set we have chosen. 

With the choice (3.14), the dPDFs involved in number sum rules but which are not 
affected by number effects satisfy their sum rules to a much better degree. It also turns out 
that once we have included terms to take account of number effects (described shortly), 
insertion of phase space factors according to (3.14) into dPDFs affected by these effects 
similarly improves the degree to which these dPDFs satisfy the sum rules. What is more, 
the momentum sum rules are much better satisfied when one uses (3.14) rather than (3.12). 
Illustration of some of these points for some representative dPDF cases, as well as an 
exposition of the extent to which we satisfy the sum rules with this choice of phase space 
factor, is given in Fig. 5. 

Having found a satisfactory phase space factor, we proceed to discuss how the second 
required feature in the list above - namely the incorporation of number effects - might be 
achieved in our input dPDFs. We have seen that number effects are particularly important 
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aT 0.6 
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Figure 5: The same sum rule ratios as in Fig. 4, but this time plotted with the phase 
space factor as in (3.14). 



for equal flavour valence-valence dPDFs, and we shall outline how suitable inputs for this 
particular type of dPDF may be constructed shortly. However, number effects can also in 
principle have an impact on any other dPDF for which the same parton type appears in 
both parton indices. Since there are only a finite number of valence up and down quarks 
in the proton (as opposed to an infinite number of sea quarks and gluons), one might 
anticipate number effects relating to these valence quarks to be most important. We now 
discuss how these effects can be included in dPDFs which 'contain' an up and/or a down 
valence combination in both of their parton indices (e.g. u + u v , d + d + , where i + = i + i). 

An example of such a distribution would be the u+u + distribution, since u + u + = 
{u v + 2u s )(u v + 2u s ), where u s = u. Consider the ways in which one can pick two up 
flavour partons (either quarks or antiquarks) from the proton. Either one can pick two sea 
partons, or one can pick a sea parton and a valence quark (in either order), or one can 
pick two valence quarks - these possibilities of course correspond to the different terms in 
the expansion of (u v + 2u s )(u v + 2u s ). Factorised terms multiplied by phase space factors 
are reasonable for all possibilities apart from the two valence option, where it would seem 
important to take account of the fact that removing a valence up halves the probability to 
find another. At a crude level we can incorporate this fact by using a term which is equal 
to half of the naive 'factorised x phase space factor' guess for the valence- valence term. We 
can think of this adjustment in another way, and say that we incorporate number effects 
in the u + u + distribution by subtracting the following term from our initial 'factorised x 
phase space factor' construct: 

^D^(x 1 -t )D^(x 2 ;t )p u ^(x 1 ,x 2 ) (3.16) 

Generalising this argument, we observe that a dPDF which contains n times the up 
valence-up valence combination in its parton indices must have n times the term (3.16) 
subtracted from it to take account of number effects. Similarly, a distribution which con- 
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tains n times the down valence-down valence combination in its parton indices must have 
n times D dv (xi;to)D d i v (x2;to)p dvdv (xi,X2) subtracted from it. Note in this case that we 
must remove the naive d v d v term entirely because there is no chance of finding two valence 
down quarks in the proton. Fig. 6 shows how inclusion of the number effect terms improves 
the extent to which dPDFs satisfy number sum rules, for a few sample cases. 

We now turn our attention to the construction of some equal flavour valence-valence 
dPDFs approximately satisfying the sum rules. The flavours we must be concerned about 
here are up, down, and strange. Note that the s v s v distribution is not zero with the given 
set of input sPDFs, even though the s v sPDF is zero. The sum rule for this dPDF reads: 

l>l — X2 

/ dx 1 D s h vSv (x 1 ,x 2 ;to) = -D 8 h + (x 2 ;to) (3.17) 
Jo 

Since the MSTW 2008LO s+ input is nonzero, the right hand side of (3.17) is nonzero, 
and consequently the s v s v dPDF cannot be zero. We can explain why the s v s v distribution 
should be nonzero by expanding the combination into double human basis pairs - s v s v = 
ss — ss — ss + ss. We expect the probability to find an ss pair to be higher than that to find 
an ss or ss pair due to number effects. Given that one has found a strange (antistrange) in 
the proton, the probability to find a further strange (antistrange) is reduced, whilst that 
to find an antistrange (strange) in addition remains the same. 

In order to construct satisfactory distributions for these three flavour types, we imagine 
that there exists a scale t < to at which only the three valence quarks in the proton may be 
resolved, and all sea distributions are zero. The sea distributions at to ar e then generated 
dynamically by DGLAP evolution between t and to- This idea has previously been put 
forward in [45-49], in which it was investigated whether the possibility exists to fit deep 
inelastic scattering data using only u v and d v inputs at a fitted low scale t. As it turns 
out, one cannot achieve a fully satisfactory fit of data using this approach, as is admitted 
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in [50]. However, since we shall only use this idea very loosely in what follows, this point 
is not of great concern to us. 

At the scale t, the only equal flavour valence-valence dPDF which can be nonzero is 
the u v u v distribution, as there is no possibility of finding two down or strange partons (be 
they quarks or antiquarks) at this scale. A suitable ansatz for the u v u v at t is a product 
of u v sPDFs multiplied by a phase space factor p appropriate at the scale, and divided by 
two to take account of valence-valence number effects: 

D^(x 1 ,x 2 ;t) = ^D^(x 1 ;t)D^(x 2 ;t)p u ^(x 1 ,x 2 ) (3.18) 

One can straightforwardly verify that the above forms for the equal flavour valence- 
valence dPDFs are consistent with the number sum rules at this scale. Now let us consider 
how the dPDFs change as we evolve from t to to under (2.1). The first two sets of terms on 
the RHS of (2.1) will mainly serve to take (3.18) into its equivalent at to (and leave the other 
equal flavour valence- valence distributions zero). However, the final set of 'sPDF feed' terms 
results in an extra contribution appearing in each equal flavour valence- valence dPDF. Only 
the —jj — jj component of an equal flavour valence-valence combination receives nonzero 
sPDF feed contributions during evolution (g — > jj contributions). Therefore, the sPDF 
feed for an equal flavour valence-valence dPDF is the following: 

_ 2 ^ (ll+I2;t) _l_P J5 (_£!_) (3,9) 

The splitting function P qg is not a very strong function of its argument (only varying 
between ^ and |). This means that, roughly speaking, we can take the sPDF feed term 
for the equal flavour valence- valence distributions as being a function of (xi + x 2 ). If we 
then ignore the subsequent effect of the first two sets of terms on the RHS of (2.1) on the 
sPDF feed contributions, then we expect the sum total sPDF feed contribution to each 
valence-valence dPDF at to to be a function of (xi + x 2 ) only: 

D{^(x u x 2 ;t ) = \~ 1 Di v (x 1 ;t )Di v (x 2 ;t )p>^(x 1 ,x 2 ) - 2g%i +a 2 ;* ) (3.20) 

We shall refer to the function g^ (x\ + x 2 ;to) as the jj correlation term, as it represents 
the 'nonf actor ised' part of the jj (or jj) distribution which is built up from correlation- 
inducing sPDF feed contributions. How should we decide on the form of this function for 
a particular choice for the flavour j? We can answer this question by using the number 
sum rule that (3.20) must satisfy, which we shall write here as: 

/•i— £2 . . . 

/ dx 1 r% 3v (x 1 ,x 2 ;to) = {N jv - l)D%(x 2 ;t ) - 2D{(x 2 -t Q ) (3.21) 
J o 

The first term on the RHS of (3.20) integrates to give approximately the first term on 
the RHS of (3.21). The —2g^{x\ + x 2 ; to) must therefore integrate to give the second term 
on the RHS of this equation: 

rl-X2 

-2/ dx 1 gii(x 1 + x 2 ;t ) = -2D] i (x 2 ;to) (3.22) 
J o 
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Figure 7: Left panel: The sum rule ratio for the u v u v number sum rule when D^ vUv is 
constructed according to (3.20) and (3.23). The ratio is close to 1 over most of the range 
of x, except near x = 0.05 where it diverges violently. This appears to indicate that the 
sum rule is being badly violated near x = 0.05. 

Right panel: The u v u v sum rule integral plotted against the sPDF quantity it should be 
equal to. This plot reveals that the divergence in the sum rule ratio is caused by the 
integral curve slightly missing a zero in the sPDF quantity, and is not serious in practice. 



This is an integral equation with a unique solution, and it is straightforward to show 
that the solution is the following: 

^ Mo) = -MgiM p.23) 

Our proposed form for the input equal flavour valence- valence distributions is therefore 
(3.20) with g^ v given by (3.23). Clearly the d v d v and s v s v number sum rules will be 
perfectly satisfied using this form. Fig. 7 shows how well the u v u v sum rule is satisfied. 

Unfortunately, with this choice for the equal flavour valence- valence dPDFs, the uu, dd, 
ss and ss dPDFs all go negative. Naively, one might view this as arising because the forms 
we have used for the equal flavour valence- valence dPDFs are in some way unsatisfactory. 
However, instead we observe that it occurs because we have omitted an important term 
in our above treatment of the distributions. Since these distributions contain the 

parton combination jj + jj that also appears in the j v j v distribution with the opposite 
sign, the receive the same sPDF feed contributions as the j v j v during evolution, but 
with the opposite sign. Thus for consistency each distribution should have an extra 
term added onto it equal to plus 2gi J (xi + X2; to). With this alteration, all double human 
basis dPDFs are again positive, and we see little adverse effect on the extent to which the 
sum rules involving j+j+ distributions are satisfied. 

Having now completed our description of how we constructed some suitable input 
dPDFs, we conclude our discussion with a short summary of how well the dPDFs satisfy 
the complete set of sum rules. In the context of the double human basis, the sum rule ratios 
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Figure 8: Left panel: The sum rule ratio for the (S + g)T% momentum sum rule, plotted 
using the fully constructed set of input dPDFs. 

Right panel: The (S+g)Ts momentum sum rule integral plotted against the sPDF quantity 
it should be equal to. 



are all within 25% of 1 for x < 0.8. Above this value, the sum rules are not obeyed so well 
- however the values of the PDFs are tiny at these x values, so large/small sum rule ratio 
values at these x values are not in practice too great a problem. In the double evolution 
basis the story is the same, barring trivial divergences due to the sum rule integral slightly 
missing a zero in the sPDF quantity it should be equal to. The one exception to this is the 
case of the T3(X + g) momentum sum rule. The sum rule ratio for this sum rule, plotted 
in the left panel of Fig. 8, plunges to 0.65 around x = 0.02. This possibly looks worse than 
it is - if one plots both the integral and the sPDF quantity it should be equal to (right 
panel of Fig. 8), then one notices that the dip in the sum rule ratio is due to the integral 
slightly overestimating a dip in the sPDF quantity in a region where the sPDF quantity is 
rather small. Furthermore, it is unlikely that the particular combination T^T. + g) will be 
directly accessed by any scattering processes at the LHC. Consequently we are prepared 
to accept the large deviation from 1 in the T^(T, + g) sum rule ratio. 

4. Numerical Solution of the Double DGLAP Equation 

There exist several options for the broad numerical method to use to integrate the dDGLAP 
equations. One could choose to adapt either the direct x space or Mellin transform methods 
which are commonly used to numerically integrate the sDGLAP equation (see, for example, 
[51,52] for routines using the x space method for solution of the sDGLAP equation, and [53] 
for a routine using the Mellin transform method). Alternatively, one could develop a 
numerical method based on the explicit solution of the dDGLAP equation in terms of 
sPDFs (2.8). This is the approach that has been preferred in the previous numerical 
treatments of the subject [37, 38]. Here we adopt an x space method. This has the 
advantages that it is conceptually simple, is flexible enough to take the inputs described in 
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Section 3.2 with no problems, and is competitive in efficiency with the other methods in 
the context of the dDGLAP equation. It also has the advantage over the 'explicit solution' 
method in that the Dj (x; t) Green's functions, which are difficult to calculate numerically 
to a sufficient degree of accuracy, do not feature. 

4.1 The dDGLAP Evolution Program 

Our program solves the dDGLAP equation (2.1) directly using a grid in xi,X2 and t. We 
choose the spacing of the grid points in t to be linear - this is the 'natural' choice, and it is 
adopted in a number of sDGLAP x-space routines (e.g. [51,52]). In the x\ and X2 directions, 
the points are taken to be evenly spaced in the variable u = ln(y^), with equal numbers of 
points in the x\ and X2 directions (600 for the grids of [40]). This gives a spacing uniform 
in ln(x) in the small x regions and directions in which the dPDF is diverging rapidly, and 
a linear spacing in larger x regions and directions in which the variation of the dPDF 
is slower. The boundary of the grid in (xi,X2) space is defined by the lines xi = x m - m , 
X2 = Xmm, xi = 1 — x m i n , X2 = 1 — x m i n , and x\ + X2 = 1 (the kinematical boundary), with 
a default x m i n = 10~ 6 . The methods we use for the numerical integration of the first two 
terms on the right hand side of the dDGLAP equations are described in the Appendix. 

The final set of terms in the dDGLAP equation (the 'sPDF feed' terms) are obtained 
at a given t by numerically evolving the sDGLAP equations contemporaneously with the 
dDGLAP equations. The grid used for the sDGLAP evolution is the similar to that used 
for the dDGLAP evolution. The only difference is that it extends in just one x direction, 
between x m i n and (1 — x m i n ). For consistency, the sPDF inputs used are the MSTW2008LO 
inputs. 

Given the structure of the dDGLAP equation, the dDGLAP evolution routine requires 
the values of the sPDFs at x values of the form Xi + Xj , where Xi and Xj are two x values 
on the uniform in ln(x/(l — x)) grid. With the grid used, it is clear that Xj + Xj does not 
also lie on the grid, so interpolation has to be used to obtain the sPDF values required. 
Away from the edges of the sPDF x-grid, natural cubic spline interpolation based on the 
sPDF values at the nearest four grid points is used, whilst linear interpolation is used at 
the edges. 

The program uses the 'double evolution' basis introduced in Section 3 as its internal 
basis for the evolution of the dPDFs. Use of this basis for the evolution is advantageous 
because the dDGLAP equations become in some sense 'minimally coupled' in this basis. 
Out of the 91 equations, 66 are rendered diagonal at LO using this basis (i.e. rate of change 
of with t is given only by the two integral terms involving D l £ , with no nonzero sPDF 
feed terms). The remaining equations have very few terms on the RHS (two terms in each 
integral term plus one sPDF feed term) . The use of this basis makes the coding in of the 
dDGLAP equations manageable. 

Stepwise evolution in t is carried out by a fourth-order Runge-Kutta method. The 
evolution begins at a scale to equal to that at which the input distributions are defined 
(Ql = 1 GeV 2 with the MSTW2008LO inputs). The final scale obtained in the evolution 
tf and the number of Runge-Kutta steps used to reach this scale Nf may be specified by 
the user. To produce the grids of [40], 120 points were used in the t direction. 
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4.2 Flavour Number Schemes 

Our program has the potential to perform the evolution using either a fixed or (zero mass) 
variable flavour number scheme, with rif fixed at 3,4,5 or 6 in the FFNS, or potentially 
varying from 3 — > 6 in the ZM-VFNS. The scheme can be determined by the user via 
the variables LGMCSQ, LGMBSQ and LGMTSQ which are equal to the thresholds in t at which 
the charm, bottom and top flavours become active respectively. For a FFNS of given nj, 
LGMCSQ, LGMBSQ and LGMTSQ should be set appropriately either above to or below tf (e.g. 
for a FFNS with n f = 5, set LGMCSQ < t ,LGMBSQ < t and LGMTSQ > tf). For a ZM-VFNS, 
at least one of LGMCSQ, LGMBSQ and LGMTSQ must lie in between to and tf. It should be 
noted that to produce the grids of [40], the program was run under a ZM-VFNS with n/ 
varying between 3 and 5. The variables LGMCSQ and LGMBSQ were set according to the 
values of m c and preferred by MSTW - 1.40 GeV and 4.75 GeV respectively. 

Prior to the evolution, the program compares LGMCSQ, LGMBSQ and LGMTSQ with to and 
tf. Depending on the results of this, it splits the full evolution from to an d tf into up to 
four intervals, each with a different value of rif. The total number of integration steps in 
t, N t , is divided up amongst these intervals roughly in proportion to the interval sizes in t. 

In each interval, the strong coupling constant t is calculated according to the LO 
analytic form: 

a ,( t ) = ^(f) . b = 33 - 2n f 

as[t) l + a s (t')b(t-t'Y 12tt ■ 1 ' 

The quantity t' corresponds to the value of t at the beginning of the interval. In the first 
interval, the boundary value of the strong coupling constant, as(t'), is taken to be the 
initial value specified by the user as (to). In later intervals it is chosen to ensure continuity 
in as, which is the appropriate matching condition at LO [41]. 

4.3 Accuracy of the Program 

We wish to get a rough estimate of the error in the dPDF values at Q introduced by 
numerical evolution with N x points in each x direction, and Nt points in the t. To do this, 
one might propose doing an evolution with twice as many points in each direction, and 
then taking the error in the original dPDFs at Q as being the absolute difference between 
the dPDF values produced by the two evolutions. Unfortunately, we cannot perform this 
procedure for the values of N x and N t used to produce the grids in [40] (600 and 120). This 
is because doubling the number of x points in this case causes the program to require far 
more RAM than a typical modern machine can provide. Instead, we show here that the 
accuracy of the program is reasonable even when N x and Nt take on the smaller values of 
150 and 10 respectively - we then know that the accuracy of the procedure with N x = 600 
and N t = 120 should be very good. 

We perform the error estimation evolution from Qo = 1 GeV to Qf = 100 GeV. In 
Fig. 9, the fractional error in the distribution Dff along the sample line x\ = X2 = x as 
calculated by the above method is plotted. That is, we plot: 

, . n v = I D h (x,x,Q f ) Nx=15 Q, Nt=10 -Dff (x,x,Qf) Nx=3 oo tNt= 2o | 

6(X,LJf) — nflfl/ r\ \ • 

D y h y (x,X,Qf) Nx=3 00,N t =20 
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Figure 9: An estimation of the numerical error when one performs an evolution from 
Q = 1 GeV to Q = 100 GeV using a grid with 150 points in each x direction, and 10 in 
the t. The error values plotted are those in the gg dPDF along the line x\ = X2 = x. 



We choose to look at Dff because this is one of the dPDFs which should be calculated 
least accurately by an evolution routine. As expected, the error increases as one approaches 
the kinematical bound due to the fact that less x points are used in the evolution inte- 
grations for the dPDF values closer to the bound. We see that the error is small in the 
crucial small x region - less than 1% for x < 0.3, and less than 6% for x < 0.4. The error 
becomes large as one approaches x = 0.5, but since this region is not likely to be important 
in applications at the LHC (which probes x\,X2 < 0.1), this is not a major problem. The 
graph indicates that even with N x = 150 and Nt = 10 the numerical evolution to LHC 
scales introduces errors which are less than 1% for x\ < 0.3, X2 < 0.3, and less than 6% for 
xi < 0.4, x 2 < 0.4. 

5. Properties of the dPDFs 

We have seen that there are two ways to improve on using simple products of sPDFs as the 
dPDFs at the (high) scale Q. First, one can use dDGLAP evolution to obtain the dPDFs 
at Q, with a reasonable choice of dPDFs at a low scale Qo used as the starting point for the 
evolution. Second, one can use improved inputs at the low scale Qo, which take account 
of momentum and number effects. In this section, we describe and illustrate the extent to 
which introducing these improvements changes the dPDFs at the scale Q. 

The large number of dPDFs precludes the possibility of discussing them all. Instead, 
we choose to focus on a small number of parton pairings which should be important in 
double scattering processes at the LHC, and which in some sense might be considered 
to form a representative set. These are the uu, uu, ug and gg pairings. Note that we 
have a dPDF for which our input form contains a valence number effect term in this set 
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Figure 10: Plots of the ratio R% input denned in equation (5.1) for Q = 100 GeV, p = 0, 1 
and 2, and the parton combinations ij discussed in the text. 



(the uu), and a distribution for which our input contains a jj correlation term (the uu). 
Furthermore, we see that the set covers all types of sPDF feed term that can appear in 
dDGLAP evolution. 

For the purposes of making concrete comparisons between different methods of obtain- 
ing the dPDFs at a high scale Q, we also need to make a specific choice for Q. Except 
where otherwise stated, we make the reasonable choice Q = 100 GeV (~ Mw,Mz, for 
example). At the scale Q, we only look at the dPDF values along the line x\ = X2 - this 
allows us to produce easily readable 2D plots. 

The main novel component of the present work is the introduction of the improved 
input dPDFs of Section 3.2. Consequently, the first question we should like to answer is how 
use of the improved inputs in the dDGLAP equation, as opposed to naive 'factorisedx (1 — 
x\ — X2) pi inputs, affects the dPDFs at the scale Q. To this end, we have plotted the 
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following ratio for our sample dPDFs in Fig. 10: 



r>ij 



(x;Q) = 



D%(x,x;Q) 



input (xi ,X2\Qo)=D" l h (x\-,Qo)D ] h (x2\Qo){l-x\-X2) p 



(5.1) 



'Ainput 




(x,x;Q) 



input 



(xi ,X2;Qo)=onv improved inputs 



We have made plots for each of the common traditional choices for p - 0, 1 and 2. 
One immediately notices in Fig. 10 that all of the ratio curves deviate significantly from 
1. This shows that the precise choice of inputs at the low scale has an important impact 
on the high scale dPDFs, and demonstrates the inadequacy of the traditional naive input 
forms. We see that multiplying factorised inputs by a phase space factor of (1 — x\ — x<i) or 
(1 — x\ — X2) 2 gives high scale dPDFs which are generally too small for small (x\, X?). This 
is expected - we have seen that (1 — x\ — X2) or (1 — x\ — X2) 2 phase space factors suppress 
the inputs too much in the high x\, low X2 and high X2, low x\ regions. Since these regions 
directly feed the small x±,X2 region, this directly translates into a deficiency in the high 
scale dPDFs in the small X\,X2 region. Conversely, we see that not using a phase factor 
in the inputs results in high scale dPDFs which are generally too large. This is because in 
this scenario the inputs are too large near the kinematic bound, and this excess propagates 
down to smaller x±,X2 values during evolution. 

It is interesting to note that, contrary to the previous general statement, the p = ratio 
for the uu dPDF actually dips below unity between x = 0.005 and x = 0.15. Furthermore, 
we see that the p = uu ratio rises above the corresponding ratios for the other flavour 
combinations. The origin of each of these features is in the extra terms we included in 
our improved inputs to take account of valence number effects or jj correlations, which 
do not appear in the naive inputs. The inclusion of a positive jj correlation term in the 
uu distribution causes our uu dPDF to be larger at the high scale than it would be if the 
correlation term were absent. Since our dPDFs appear on the denominator of R l ^ input , 
this manifests itself as a reduction in our p = uu ratio. Conversely, the subtraction of a 
valence number effect term from our uu input results in a reduction of our uu dPDF at Q, 
which increases the uu ratio. 

For p = 1 and 2, we observe that the uu ratio is still larger than the others for small 
x. However, the uu ratio is now very slightly larger than the ug and gg ratios at small x 
values. This is because the ug and gg high scale distributions at small x are more sensitive 
to the form of the input distributions near the kinematic boundary than the uu. This is a 
simple consequence of the fact that gluon type evolution causes a faster cascade of PDFs 
to low x values than u or u type evolution. The reduction in the ug and gg ratios at small 
x relative to the uu due to the change in p overcomes the small effect of including the jj 
correlation term in our uu. 

The contributions of the jj correlation and valence number effect terms to the high 
scale (Q = 100 GeV) double human basis dPDFs are most cleanly observed at x ~ 0.05, 
and are on the order of 10% in this x region. For smaller x, the contributions from the 
extra terms are swamped by sea-sea contributions to the dPDF, whilst at larger x, phase 
space effects become dominant. 



- 29 - 



c 
S 
o 

ft? 



1.2 
1 

0.8 
0.6 
0.4 
0.2 




uu 



Q = IGeV 
Q = lOGeV 
Q = lOOGeV 
Q = lOOOGeV 



10~ 6 10~ 5 IO" 4 10~ 3 10~ 2 IO" 1 10° 

x 



ft? 



1.2 
1 

0.8 

0.6 \ 

0.4 

0.2 


10 



— i . — n — ■ ■ — n — ■ ■— n — ■ ■ — n — ■ n — ■ ■ — *■ 

ug 




\ s \ 
\ x \ 






Q = IGeV 

. Q = lOGeV 

Q = lOOGeV 

Q = lOOOGeV 

- - ■ - - ■ - - ■ 





io~ 5 io- 4 io-3 io~ 2 io- 1 10° 



ft? 



1.2 
1 

0.8 
0.6 
0.4 
0.2 




Q = IGeV 
Q = lOGeV 
Q = lOOGeV 
Q = lOOOGeV 



10~ 6 10~ 5 IO" 4 10~ 3 10~ 2 IO" 1 10° 

x 



ft? 



1.2 
1 

0.8 
0.6 
0.4 
0.2 




— i . — n — ■ ■— n — i ■ i — ■ ■— n — ■ ■ — n — ■ ■ — *■ 

99 










\ v \ 




\ \ 
\ V \ 

'■ * \ 
\ 1 \ 
' \ \ 


Q = IGeV 

. Q = lOGeV 

Q = lOOGeV 

Q = lOOOGeV 

- - ■ - - ■ - - ■ 





10 



>- 6 io~ 5 io- 4 io~ 3 io~ 2 io- 1 10° 



Fi gure 11: Plots of the ratio R^ yo £ CC( j defined in equation (5.2) for Q — 1, 10, 100 and 
1000 GeV and the parton combinations ij discussed in the text. 



Aside from looking at the effect of using different inputs on the dPDFs at scale Q, we 
can also ask to what extent correlations introduced by dDGLAP evolution affect the dPDFs 
at Q. There are essentially two types of correlations that the dDGLAP equation introduces 
- correlations due to the requirement of momentum conservation, and more interesting 
correlations generated by the sPDF feed terms. Here, we choose to look specifically at the 
effect of the latter. 

In order to do this, we evolved our improved input dPDFs up to the scales Q = 10 GeV, 
Q = 100 GeV, and Q = 1000 GeV, both with the sPDF feed terms included in the evolution, 
and also with these terms set to zero. For each final scale and parton pairing in our selected 
set, the following ratio was then plotted: 



T-.ii / (x, X, \our improved inputs, no sPDF feed i~ ~.\ 

R L feedfo Q) = ~iTr ^Ti ( 5 - 2 ) 



1)^ (x, X] Q) | our improved inputs 



- 30 - 



We plot the results using a logarithmic x scale in Fig. 11 . The effect of the sPDF 
terms is small but non-negligible, being at roughly the 10% level for x < 10 -2 in all of the 
dPDFs considered, and increasing with Q. 

We observe that the ratios for all of the given flavour combinations look very similar 
for x from 1CT 6 to 10~ 4 . The reason for this is that the small x shape of the distributions 
considered is very strongly determined by the (either direct or indirect) feeding of these 
distributions by the gg distribution. If the gg dPDF loses its sPDF feed and is reduced by 
a certain percentage at small x, the connection of the other dPDFs to the gg will result 
in these dPDFs being reduced by a similar amount. This explanation can be verified by 
investigating what happens if we remove all of the sPDF terms except for the gg feed. In 
this case the ratios for all of the considered dPDFs are much closer to 1 for 10 -6 < x < 10~ 4 , 
suggesting that the subtraction of the gg sPDF feed is the dominant factor determining 
the shapes of the plots in Fig. 11 for small x. 

For larger x, the deviation of the uu ratio from 1 remains small, and tends to as x 
approaches its maximum of 0.5. This is expected since there is no direct sPDF feed term 
in the evolution of the uu dPDF. The uu ratio also seems to tend to 1 as x — > 0.5, albeit 
more slowly, whilst the ug and gg ratios plunge towards zero, the gg more rapidly than the 
ug. This implies that at large x, the sPDF feed contributions are more important to the gg 
than they are to the ug, and that they are more important to the ug than they are to the 
uu. We can explain this ordering using a fact we have previously mentioned - namely, that 
the 'pull' on a gluon PDF towards lower x values during evolution is stronger than that on 
a quark type PDF. The gg distribution at large x is pulled strongly towards lower x values 
in two directions, and is very much smaller if it is not continuously fed by an sPDF. By 
contrast, the 'pull' on the large x uu distribution is smaller in both directions, and so the 
contribution of similar sPDF feed terms is proportionately smaller. The ug distribution 
has one gluon flavour index and one quark, so the importance of the sPDF feed on this 
distribution at large x is intermediate. 

We have not been able to exactly reproduce the results of either of the extant numerical 
investigations into the correlations induced by evolution - [37] and [38]. However, we do 
agree with [38] that the accumulated sPDF feed contribution to the gg between ~ 1 GeV 
and 100 GeV accounts for about 10% of the Q = 100 GeV gg distribution at small x. In 
Fig. 12, we plot the following ratio for Q = 80.4 GeV: 

factoriscd inputs 

-Dj{x-Q)D{{x-Q) 
( ' Q)= Di(x- 1 Q)Di(x;Q) (5 " 3) 

This figure corresponds to the solid curve in Fig. 1 of [37], with MSTW2008LO inputs 
replacing the MRS99 inputs used there. We expect that the ratio R" should tend to — 1 
as x approaches 0.5 for any Q sufficiently larger than the input scale. This is because 
evolution will very quickly cause Dff to become much smaller than the factorised value 
near the kinematic bound. Our curve exhibits this property, but it seems unlikely that the 



7 In this figure, and in figures 12 and 14, we make plots down to x = 1CP 6 . Although it is interesting 
to look at our LO dPDFs at very small x, we should mention that we do not expect the leading order 
approximation to produce very accurate dPDFs in this region. 
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Figure 12: gg correlation ratio B?o at Q = 80.4 GeV obtained using MSTW2008LO 
f actor ised inputs. 



solid curve plotted in Fig. 1 of [37] will, especially if it reaches 0.6 for higher x values as is 
stated in [37]. 

Finally, we compare our full treatment (improved inputs plus full dDGLAP evolution) 
with the approximation that simply uses factorised inputs x(l — x\ — X2Y (p = 0, 1 or 
2) at the scale Q. This approximation is frequently used in phenomenological studies of 
double parton scattering processes. In Fig. 13, we plot the following ratio along the line 
x\ = X2 = x for our sample dPDFs and for p = 0, 1 and 2: 



D i h (x 1 ;Q)Di(x 2 ;Q)(l - Xl - x 2 )p 
D%(x 1 ,x 2 ;Q) I our improved inputs 



(5.4) 



The plots reveal that even a (1 — x\ — x 2 ) 2 phase space factor multiplying a factorised 
form at Q underestimates the large x falloff in the dPDFs along x\ = x 2 = x. For very 
small x, the ratios are all slightly less than 1 due to the fact that one misses the sPDF 
feed contributions if one uses a factorised form at Q (note that the ratio appears smallest 
at very low x for the uu, due to the fact that the sPDF feed for the uu is particularly 
important around x = 10 -2 - see Fig. 11). One also notices the imprint of omitting the 
valence number effect and jj correlation terms in the ratios - the uu ratio rises above the 
others at x ~ 0.05, whilst the uu dips at this x value. 

It is interesting to consider the behaviour of R l j\f ina i{xi-,x 2 ;Q) away from the line 
x\ = x 2 = x. In Fig. 14, we plot the p = ratio for the gg flavour combination along 
several lines emanating from the point x\ = 10~®,x 2 = 10~ 6 . The figure shows that the 
deviation of this ratio from 1 is maximal along x\ = x 2 (in fact, this statement holds for 
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Figure 13: Plots of the ratio ^X/inoi defined i n equation (5.4) at Q = 100 GeV and along 
the line x\ = X2 = x. The ratio is plotted for p = 0, 1 and 2 and for each of the parton 
combinations ij discussed in the text. 



any combination of parton indices). We observe that a p = factorised form is a fairly 
good approximation to our gg dPDF close to the x\ axis, except when x\ is very large 
(#i > 0.8). This is to be expected, given our use of input dPDFs which essentially reduce 
to p = factorised forms near the lines x\ = and X2 = 0. One can infer from the plot 
that use of a factorised form multiplied by either (1 — x\ — x%) or (1 — x\ — X2) 2 will result 
in one overestimating the falloff in the dPDFs in the x\ ~ 0,a?2 < 0.8 and x\ < 0.8, X2 ~ 
regions. 

6. Summary and Outlook 

In this report, we have developed a framework based on the dDGLAP equation for cal- 
culating the LO double distributions DV(x±,X2',t) which represents an improvement on 
approaches used previously. We have derived for the first time the momentum and sum 
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plotted along various lines of the form X2 
100 GeV. 



rules that the dPDFs have to obey. An important implication of these sum rules is that the 
conventionally held wisdom that the dPDFs should be approximately equal to products 
of sPDFs for small x\,X2 does not apply in the case of the equal flavour valence- valence 
dPDFs. Using the dPDF sum rules, we have constructed a set of improved input dPDFs 
corresponding to the MSTW2008LO sPDF inputs. In the double human flavour basis, 
these dPDFs are all positive and satisfy the sum rules to better than 25% precision (for 
x < 0.8). 

We have written a program which numerically integrates the LO dDGLAP equation 
using a direct x space method, enabling one to evolve the dPDF inputs to higher scales. The 
accuracy of the program is good for small xi,X2 - an evolution from 1 GeV to 100 GeV 
using a grid with only 150 points in each x direction and 10 points in the t direction 
produces dPDF values with numerical errors of less than 1% for x\ < 0.3, X2 < 0.3. We 
have produced a set of publicly available dPDF grids by applying the numerical procedure 
to our improved inputs, which can be found along with interpolation code at [40]. To 
produce the grids, 600 points were used in each x direction, and 120 in the t, ensuring an 
accuracy much better than 1% for small x. 

We saw that the accuracy of our program is rather poor near the kinematical bound. If 
the accuracy here needed to be improved without significantly increasing computing time, 
then a multigrid method could be implemented in the program (for an example of the use 
of this method for the sPDF case, see [51]). The additional more finely spaced grids would 
be introduced in the region near the kinematic bound to increase accuracy in this region. 

For the purposes of experimental studies of double parton interactions in the near 
future, which will be attempting to establish the existence of correlations in the dPDFs, 
the LO treatment presented here is sufficiently accurate. If correlations are found and they 
agree with some or all of the predictions made here, then this will be a strong impetus for us 



to extend the formalism to NLO. As mentioned in Section 2, such an extension is not trivial, 
since the structure of the third term of the dDGLAP equation becomes significantly more 
complex at NLO. At this order, one requires the functions P^j k (x±,X2) which cannot be 
obtained trivially from the NLO sPDFs as in the LO case. It is likely that these functions 
exist in the literature, although some work may need to be done to get them into a form 
that can be used in the dDGLAP equation. 

Many double scattering processes which might provide important signals/backgrounds 
at the LHC do not involve the same hard scale in both collisions. An example of such a 
process is the simultaneous production of a W and a bb pair in separate collisions. This 
forms a background to the process p + p — > WH, H — > bb, which might be an important 
process to discover the Higgs if < 2mw [20]. To make theoretical predictions relating 
to these processes, we require the more general double distributions with t± 7^ t<i- As 
is mentioned in Section 2, we believe that these distributions are calculated by adding 
an extra sDGLAP evolution in one variable on top of the dDGLAP evolution. A useful 
extension to the work would be to produce a more general set of double distributions based 
on this hypothesis. 

Finally, there exists the possibility of using the dPDFs developed above to undertake 
a phenomenological investigation of double parton scattering at the LHC. In particular 
it would be interesting to examine how the 'correlations' introduced via our inputs and 
by evolution affect the properties of a double scattering event, and also how one might 
measure the correlations in practice. We are currently in the process of making such an 
investigation. 
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Appendix 

A. Numerical techniques for evaluating the dDGLAP integrals 

Let us consider the integrals which have to be numerically approximated using the (xi, x^) 
grid. All of these integrals are of the following schematic form: 

I(y)= I' " -D(z,y)P (-) (A.l) 

Jx z KzJ 

The splitting function P(x) may in general consist of three terms. The first of these is 
a regular term A(x) and the second is a term proportional to a delta function KS(1 — x). 
The final term consists of a product of two factors. The first of these is a simple regular 
function R(x), whilst the second is a function S(x) containing a singular factor 1/(1 — x) 
which is regularised by the plus prescription: 

P{x) = A(x) + K5(l -x) + R(x)[S(x)]+. (A.2) 

Inserting the form (A.2) into (A.l), we find that the integrals which have to be approxi- 
mated using the grid have the following general form: 
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= h(y) + h{y) + h{y) with 
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(A.6) 



The integral in the last term of (A.6) can be done analytically for each splitting function. 
The integrals in (A.4) and the first term of (A.6) are the ones that must be performed 
on the grid. We note that the integrand in the first term of (A.6) has the property that 
it is undefined for z = x (due to the fact that S(x/z) contains a factor 1/(1 — x/z)). It 
nevertheless tends to a finite limit as z — > x (due to the fact that the divergence in S(x/z) is 
compensated for by the other factor in the integrand going to zero as z — > x). This suggests 
the use of a method for performing the numerical integrations which effectively estimates 
the integrand between z = x and the grid point with next highest z by extrapolating from 
integrand values on nearby grid points (with higher z). 

A method which uses an open Newton-Cotes rule of degree n for the first n integration 
intervals, and then switches to a closed Newton-Cotes rule to perform the integration over 
the remaining intervals, has this property. If the number of integration intervals is greater 
than 3, we use Simpson's rule as the closed rule, combined with an open rule of degree 4 
when the number of integration intervals is even, and an open rule of degree 5 otherwise. 
Open rules of the appropriate degree are used on their own when the number of intervals 
is 3 or fewer. This ensures an overall integration method which for most integrals has an 
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error of O (n An 5 d ^JiP ) ■ In this formula, n is the number of intervals used, An is the 
(even) grid spacing in u = ln(^-§-), / is the integrand taking into account the Jacobian on 
the transformation into u space, and £ is the value of u that maximises df^/d 



u 



With the numerical method described, the integral (A.l) is approximated by: 
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(A.7) 



The indices {i,j,k} represent grid points, with i corresponding to the grid point with z 
value equal to x {z{ = x) and k corresponding to the point with z value equal to 1 — y 
{zk = 1 — y). The w^k are Newton-Cotes type integration weights whose values are 
dictated by the prescription described above. Note that the weight at grid point j under 
this prescription depends on the start and end points of the integration - hence w depends 
on the indices i and k. The function J(x) is the Jacobian, J(x) = dx/du = x(l — x). 
We may rewrite (A.7) as: 



where 
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otherwise. 



(A.9) 



The three-dimensional array P^k only depends on the splitting function P(x), Jacobian 
J{x) and weights w^k- None of these vary during an evolution, with the possible exception 
of P gg (this contains a term proportional to nj in the K5(l — x) piece and so may vary 
in a variable flavour number scheme - see Section 4.2). We therefore precalculate and 
store the elements of Pijk during program initialisation, to increase efficiency The possible 
variation of the contributions to P^k from the term in P gg proportional to nj is handled 
by postponing the calculation of these contributions such that they are calculated and 
reintroduced at each evolution step (using the value of n f appropriate to that step) . 
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